\documentclass{article}
\usepackage{amsmath,amsfonts,geometry,verbatim,graphicx}
\geometry{margin=1in}
\renewcommand{\tilde}{\widetilde}
\newcommand{\Lie}{\mathcal{L}}
\newcommand{\ura}{\mathbf}
\newcommand{\includedraw}[1]{\includegraphics[width=\linewidth]{#1.drawit/QuickLook/Preview.jpg}}
\newcommand{\BLambda}{{\sf \Lambda\hspace*{-0.9ex}%
\rule{0.15ex}{1.5ex}\hspace*{0.9ex}}}
\title{Boundaries and the Conservation of Energy}
\author{Michael Barriault}

\begin{document}

\section{Source case}

Consider one of the simplest one-dimensional time/space wave equation
\begin{gather}
\dot u = -\nabla u \label{swave}
\end{gather}
where $\dot~ \equiv \partial_t$ and $\nabla \equiv \partial_s$, in the region $s \in \left(0,1\right)$ and $t \geq 0$. We will assume that $u(0,s)$ is known. The solution $u(t,s)$ comprises a wave that moves forward in time and to the right in space. Figure \ref{wavepath} shows the path travelled through spacetime by a wave crest. Note that those crests on the left are influenced by the boundary $u(t,0)$. However, $u(t,0)$ is a point like any other, and so is influenced by points {\em outside} the boundary. However, the region is an isolated box, and since waves carry energy, if a wave is affected by points outside the box it must be carrying energy from outside the box to the inside - a violation of conservation of energy. As a fix, we can say that the energy along $u(t,0)$ is at a constant value across the boundary - that is, $\left.\nabla E\right|_{s=0} = 0$. Or, since energy can be expressed with respect to some reference value, we can say without loss of generality that $\left.E\right|_{s=0} = 0$. Since energy is proportional to speed ($E \propto \dot u^2$), this is equivalent to setting $\dot u(t,0) = 0$. 

For left-travelling waves ($\dot u = \nabla u$) the same treatment is applied to the right boundary, so $\dot u(t,1)=0$.

\begin{figure}[h!]
\centering
\includedraw{wavepath}
\caption{Contours of wave crests as they travel through time and space.}
\label{wavepath}
\end{figure}

\section{System of equations}

For more complicated wave equations, such as
\begin{gather}
\ddot u = \nabla^2 u
\end{gather}
we must do a little more work. First, we must convert this to a system of first-order equations. Do this by setting $f = \dot u$ and $g = \nabla u$. Then the second-order equation becomes
\begin{gather}
\dot u = v \\
\dot f = \ddot u = \nabla^2 u = \nabla g \\
\dot g = \nabla \dot u = \nabla f
\end{gather}
We can rewrite this into matrix form
\begin{gather}
\left( \begin{array}{c} \dot u \\ \dot f \\ \dot g \end{array} \right) 
= \left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{array} \right) 
\left( \begin{array}{c} \nabla u \\ \nabla f \\ \nabla g \end{array} \right) 
+ \left( \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right) 
\left( \begin{array}{c} u \\ f \\ g \end{array} \right) 
\end{gather}
or
\begin{gather}
\dot{\ura u} = \mathbb{S}\nabla \ura u + \mathbb{O} \ura u
\end{gather}
where
\begin{gather}
\ura u = \left( \begin{array}{c} u \\ f \\ g \end{array} \right) \\
\mathbb{S} = \left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{array} \right) \\
\mathbb{O} = \left( \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right) 
\end{gather}
We wish to transform $u,f,g$ into three new functions - linear combinations of $u,f,g$ - such that each one is a simple wave described by \eqref{swave}, at least on the boundaries $s=0$ and $s=1$ where the $s$ derivative - requiring points outside the boundary - might draw in extra energy. We know from linear algebra that, if $\mathbb S$ is diagonalizable, then there is a diagonal matrix $\mathbb L$ comprised of the eigenvalues of $\mathbb S$, and an orthogonal matrix $\mathbb P$ with columns comprised of the eigenvectors of $\mathbb S$, in the same order as the associated eigenvalues in $\mathbb L$, so that
\begin{gather}
\mathbb L = \mathbb P^{-1} \mathbb S \mathbb P
\end{gather}
or
\begin{gather}
\mathbb S = \mathbb P \mathbb L \mathbb P^{-1}
\end{gather}
then if $\dot{\ura u} = \mathbb S \nabla \ura u$ (we don't care about the non-derivative terms), setting $\ura v = \mathbb P^{-1} \ura u$, then
\begin{align}
\mathbb P^{-1} \dot{\ura u} &= \mathbb P^{-1} \mathbb S \nabla \ura u \\
\dot{\left( \mathbb P^{-1} \ura u \right)} - \dot{\mathbb P}^{-1} \ura u &= \mathbb P^{-1} \mathbb P \mathbb L \mathbb P^{-1} \nabla \ura u \\
\dot{\ura v} &= \mathbb L \left( \nabla\left(\mathbb P^{-1} \ura u\right) - \nabla \mathbb P^{-1} \ura u \right) + \dot{\mathbb P}^{-1} \ura u \\
\dot{\ura v} &= \mathbb L \nabla \ura v - \left( \mathbb L \nabla \mathbb P^{-1} - \dot{\mathbb P}^{-1} \right) \ura u \\
\dot{\ura v} &= \mathbb L \nabla \ura v - \left( \mathbb L \nabla \mathbb P^{-1} - \dot{\mathbb P}^{-1} \right) \mathbb P \ura v \\
\dot{\ura v} &= \mathbb L \nabla \ura v + \left( \mathbb L \mathbb P^{-1} \nabla \mathbb P - \mathbb P^{-1} \dot{\mathbb P} \right) \ura v \\
\end{align}
The collective $\ura v$ terms on the right-hand side are fluff that involve derivatives of the already known quantity $\mathbb P$, that we merely simplified to look nice. What's important is that our equation has effectively been transformed into
\begin{gather}
\dot{\ura v} = \mathbb L \nabla \ura v
\end{gather}
which is almost identical to what we started with, except $\mathbb L$ is diagonal. If we write this in component form,
\begin{gather}
\dot v_i = L_{ij} \nabla v_j = L_{ii} \nabla v_i = \lambda_i \nabla v_i
\end{gather}
which shows that each $v_i$ has a form exactly like \eqref{swave}. The eigenvalues of $\mathbb S$ are solutions to
\begin{gather}
-\lambda\left(\lambda^2-1\right) = 0 \Rightarrow \lambda = 0, 1, -1
\end{gather}
The associated eigenvectors are
\begin{gather}
\ura e_0 = \left( \begin{array}{c} 1 \\ 0 \\ 0 \end{array} \right) \\
\ura e_1 = \left( \begin{array}{c} 0 \\ 1 \\ 1 \end{array} \right) \\
\ura e_{-1} = \left( \begin{array}{c} 0 \\ -1 \\ 1 \end{array} \right)
\end{gather}
so that the diagonalization matrix is
\begin{gather}
\mathbb P = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & -1 \\ 0 & 1 & 1 \end{array} \right)
\end{gather}
Now we can consider the functions
\begin{gather}
\ura v = \left( \begin{array}{c} v \\ u^+ \\ u^- \end{array} \right) = \mathbb P^{-1} \ura u = \left( \begin{array}{c} u \\ \frac12 \left( f + g \right) \\ \frac12 \left( g - f \right) \end{array} \right)
\end{gather}
Note that
\begin{gather}
\dot u^+ = \frac12 \left( \dot f + \dot g \right) = \frac12 \left( \nabla g + \nabla f \right) = \nabla u^+ \\
\dot u^- = \frac12 \left( \dot g - \dot f \right) = \frac12 \left( \nabla f - \nabla g \right) = -\nabla u^-
\end{gather}
So the functions $u^+$ and $u^-$ are described exactly like \eqref{swave}, where $u^+$ is travelling to the left and $u^-$ is travelling to the right. Thus we will set $\dot u^-(t,0) = 0$ and $\dot u^+(t,1)=0$. But note that
\begin{gather}
f = u^+ - u^- \\
g = u^+ + u^-
\end{gather}
so
\begin{gather}
\dot f(t,0) = \dot u^+(t,0) - \dot u^-(t,0) = \dot u^+(t,0) = \frac12 \left( \nabla g(t,0) + \nabla f(t,0) \right) \\
\dot f(t,1) = \dot u^+(t,1) - \dot u^-(t,1) = -\dot u^-(t,1) = \frac12 \left( \nabla g(t,1) - \nabla f(t,1) \right) \\
\dot g(t,0) = \dot u^+(t,0) + \dot u^-(t,0) = \dot u^+(t,0) = \frac12 \left( \nabla g(t,0) + \nabla f(t,0) \right) \\
\dot g(t,1) = \dot u^+(t,1) + \dot u^-(t,1) = \dot u^-(t,1) = \frac12 \left( \nabla f(t,1) - \nabla g(t,1) \right)
\end{gather}
Note that on the left-hand side is the values we are manually setting, whereas the right-hand side involves values with are calculating with traditional methods.

\section{Two-dimensional systems}

Now consider the two-dimensional system
\begin{gather}
\ddot u = \nabla^2 u
\end{gather}
where $u = u(t,x,y)$. We break this down as before into a system of coupled one-dimensional equations
\begin{gather}
\dot u = f \\
\dot f = g_x + h_y \\
\dot g = f_x \\
\dot h = f_y
\end{gather}
or
\begin{gather}
\dot{\ura u} = \mathbb{X}\ura u_x + \mathbb{Y}\ura u_y + \mathbb{O}\ura u
\end{gather}
The non-trivial eigenvectors (those with non-zero eigenvalues) associated with $\mathbb X$ are
\begin{gather}
\ura e_1 = \left( 0,1,1,0 \right) \\
\ura e_{-1} = \left( 0,1,-1,0 \right)
\end{gather}
which results in the linearly combined functions
\begin{gather}
U = f + g \\
V = g - f
\end{gather}
with assoicated equations
\begin{gather}
\dot U = U_x \\
\dot V = -V_x
\end{gather}
or
\begin{gather}
f = U-V \\
g = U+V
\end{gather}
Since we want to set $\dot U = 0$ where $x=1$ and $\dot V = 0$ where $x=0$, we must apply the kludge
\begin{gather}
\dot f = \dot U = \frac12 \left( f_x + g_x + h_y \right) \\
\dot g = \dot U = \frac12 \left( f_x + g_x + h_y \right)
\end{gather}
at $x=0$ and
\begin{gather}
\dot f = \frac12 \left( -f_x + g_x + h_y \right) \\
\dot g = -\frac12 \left( g_x - g_x - h_y \right)
\end{gather}
along $x=1$. The exact same treatment can be applied to $\mathbb Y$ and the boundaries $y=0$ and $y=1$, which will result in
\begin{gather}
y=0 ~ \begin{cases} \dot f = \frac12 \left( f_y + g_x + h_y \right) \\ \dot h = \frac12 \left( f_y + g_x + h_y \right) \end{cases} \\
y=1 ~ \begin{cases} \dot f = \frac12 \left( -f_y + g_x + h_y \right) \\ \dot h = \frac12 \left( f_y - g_x - h_y \right) \end{cases} 
\end{gather}

However, what do we do about the corners $(0,0),(0,1),(1,0),(1,1)$? We need to back up a few steps. Above, what was causing trouble was a wave travelling in a direction normal to the boundary coming from outside to inside. At the corners, there isn't a well-defined normal (since the sharp turn would make the normal discontinuous). But we can define a ``good enough'' normal as the sum (really the average but we don't care about the size of $\ura n$) of the limits of the normals as they approach the corner. For example, at $(0,0)$,
\begin{gather}
\ura n = \lim_{y\to 0} \left.\ura n\right|_{x=0} + \lim_{x\to 0} \left.\ura n\right|_{y=0} = \left(1,0\right) + \left(0,1\right) = \left(1,1\right)
\end{gather}
Since we want to consider derivatives in the direction of the normal, we should change to a new coordinate system such that the normals point in the direction of a coordinate (given the rectangular nature of our region, the other normals will either be $90^\circ$ to the one at $(0,0)$ or antiparallel to it). The transformation
\begin{gather}
x = \frac1{\sqrt2} \left( r-s \right) \\
y = \frac1{\sqrt2} \left( r+s \right)
\end{gather}
is satisfactory, while being orthonormal and $(r,s)$ form a right-handed coordinate system. Clearly, (in $xy$-coordinates) $\nabla r = \frac1{\sqrt2} (1,1)$ (the normal at $(0,0)$) and $\nabla s = \frac1{\sqrt2} (-1,1)$ (the normal at $(0,1)$). In these new coordinates,
\begin{gather}
\dot u = f \\
\dot f = g_x + h_y = \frac1{\sqrt2} \left( g_r - g_s + h_r + h_s \right) \\
\dot g = f_x = \frac1{\sqrt2} \left( f_r - f_s \right) \\
\dot h = f_y = \frac1{\sqrt2} \left( f_r + f_s \right)
\end{gather}
However, for convenience, we'll do another function change
\begin{gather}
g = \frac1{\sqrt2} \left( j - k \right) \\
h = \frac1{\sqrt2} \left( j + k \right)
\end{gather}
so that
\begin{gather}
\dot u = f \\
\dot f = j_r + k_s \\
\dot j = f_r \\
\dot k = f_s
\end{gather}
or
\begin{gather}
\dot{\ura u} = \mathbb{R}\ura u_r + \mathbb{S}\ura u_s + \mathbb{O}\ura u
\end{gather}
which has the exact same form as original system, except with functions $u,f,j,k$ of $r,s$ instead of $u,f,g,h$ of $x,y$. Thus, we already know how we impose energy conservation at the troublesome boundaries, in this case at the points $(0,0)$ and $(\sqrt2,0)$ for $r$-derivatives, and $\left(\frac1{\sqrt2},-\frac1{\sqrt2}\right)$ and $\left(\frac1{\sqrt2},\frac1{\sqrt2}\right)$ for $s$-derivatives.
\begin{gather}
r=0 ~ \begin{cases} \dot f = \frac12 \left( f_r + j_r + k_s \right) \\ \dot j = \frac12 \left( f_r + j_r + k_s \right) \end{cases} \\
r=\sqrt2 ~ \begin{cases} \dot f = \frac12 \left( -f_r + j_r + k_s \right) \\ \dot j = \frac12 \left( f_r - j_r - k_s \right) \end{cases} \\
s=-\frac1{\sqrt2} ~ \begin{cases} \dot f = \frac12 \left( f_s + j_r + k_s \right) \\ \dot k = \frac12 \left( f_s + j_r + k_s \right) \end{cases} \\
s=\frac1{\sqrt2} ~ \begin{cases} \dot f = \frac12 \left( -f_s + j_r + k_s \right) \\ \dot k = \frac12 \left( f_s - j_r - k_s \right) \end{cases} 
\end{gather}
or back in our original $(x,y)$ coordinates,
\begin{gather}
(0,0) ~ \begin{cases} \dot f = \frac12 \left( \frac1{\sqrt2} \left( f_x + f_y \right) g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x - f_y + \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( -f_x + 3f_y + \sqrt2 \left( g_x + h_y \right) \right) \end{cases} \\
(1,1) ~ \begin{cases} \dot f = \frac12 \left( -\frac1{\sqrt2} \left( f_x + f_y \right) + g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x - f_y - \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( -f_x + 3f_y - \sqrt2 \left( g_x + h_y \right) \right) \end{cases} \\
(0,1) ~ \begin{cases} \dot f = \frac12 \left( \frac1{\sqrt2} \left( f_y - f_x \right) + g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x + f_y - \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( f_x + 3f_y + \sqrt2 \left( g_x + h_y \right) \right) \end{cases} \\
(1,0) ~ \begin{cases} \dot f = \frac12 \left( \frac1{\sqrt2} \left( f_x - f_y \right) + g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x + f_y + \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( f_x + 3f_y - \sqrt2 \left( g_x + h_y \right) \right) \end{cases}
\end{gather}

The following table geometrically demonstrates all the required energy conservation fixes along the boundary of the region.

\begin{tabular}{r|c|l}
$\begin{cases} \dot f = \frac12 \left( \frac1{\sqrt2} \left( f_y - f_x \right) + g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x + f_y - \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( f_x + 3f_y + \sqrt2 \left( g_x + h_y \right) \right) \end{cases}$ & 
$\begin{cases} \dot f = \frac12 \left( -f_y + g_x + h_y \right) \\ \dot g = f_x \\ \dot h = \frac12 \left( f_y - g_x - h_y \right) \end{cases} $ & 
$\begin{cases} \dot f = \frac12 \left( -\frac1{\sqrt2} \left( f_x + f_y \right) + g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x - f_y - \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( -f_x + 3f_y - \sqrt2 \left( g_x + h_y \right) \right) \end{cases}$ \\ \hline
$\begin{cases}
\dot f = \frac12 \left( f_x + g_x + h_y \right) \\
\dot g = \frac12 \left( f_x + g_x + h_y \right) \\
\dot h = f_y
\end{cases}$ & 
&
$\begin{cases}
\dot f = \frac12 \left( -f_x + g_x + h_y \right) \\
\dot g = -\frac12 \left( g_x - g_x - h_y \right) \\
\dot h = f_y
\end{cases}$ \\ \hline 
$\begin{cases} \dot f = \frac12 \left( \frac1{\sqrt2} \left( f_x + f_y \right) g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x - f_y + \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( -f_x + 3f_y + \sqrt2 \left( g_x + h_y \right) \right) \end{cases}$ &
$\begin{cases} \dot f = \frac12 \left( f_y + g_x + h_y \right) \\ \dot g = f_x \\ \dot h = \frac12 \left( f_y + g_x + h_y \right) \end{cases}$ &
$\begin{cases} \dot f = \frac12 \left( \frac1{\sqrt2} \left( f_x - f_y \right) + g_x + h_y \right) \\
\dot g = \frac14 \left( 3f_x + f_y + \sqrt2 \left( g_x + h_y \right) \right) \\
\dot h = \frac14 \left( f_x + 3f_y - \sqrt2 \left( g_x + h_y \right) \right) \end{cases}$
\end{tabular}

\section{Maxwell's equations}

We will now apply this procedure to Maxwell's (homogeneous) equations inside the unit cube
\begin{gather}
\dot{\ura E} = \nabla \times \ura B \\
\dot{\ura B} = -\nabla \times \ura E
\end{gather}
Here we have 6 coupled equations, six boundaries that involve derivatives in one coordinate direction (the faces), twelve boundaries that involve derivatives in two coordinate directions (the edges), and eight boundaries that involve derivatives in three coordinate directions (the vertices).

Let
\begin{gather}
\ura F = \left( \begin{array}{c} E^x \\ E^y \\ E^z \\ B^x \\ B^y \\B^z \end{array} \right) \\
\mathbb X = \left( \begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 0 \end{array} \right) \\
\mathbb Y = \left( \begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \end{array} \right) \\
\mathbb Z = \left( \begin{array}{cccccc} 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{array} \right)
\end{gather}
then
\begin{gather}
\dot{\ura F} = X \ura F_x + Y \ura F_y + Z \ura F_z
\end{gather}
The eigenvectors of $\mathbb X$ give four functions of interest
\begin{gather}
G = \frac12 \left( B^y + E^z \right) \\
H = \frac12 \left( B^z - E^y \right) \\
J = \frac12 \left( B^z + E^y \right) \\
K = \frac12 \left( B^y - E^z \right)
\end{gather}
where (denoting $\sim$ as terms involving other derivatives which we don't particularly care about at the moment)
\begin{gather}
\dot G = \frac12 \left( \dot B^y + \dot E^z \right) = \frac12 \left( E^z_x - E^x_z + B^y_x - B^x_y \right) = G_x + \sim \\
\dot H = H_x + \sim \\
\dot J = -J_x + \sim \\
\dot K = -K_x + \sim 
\end{gather}
Thus we want to set $\dot J = 0$ and $\dot K = 0$ on the surface $x=0$, and $\dot G = 0$ and $\dot H = 0$ on the boundary $x=1$. Rearranging,
\begin{gather}
E^y = J-H \\
E^z = G-K \\
B^y = G+K \\
B^z = J+H
\end{gather}
so now
\begin{gather}
x = 0 ~ \begin{cases}
\dot E^y = \frac12 \left( E^y_x - B^z_x - E^x_y + B^x_z \right) \\
\dot E^z = \frac12 \left( E^z_x + B^y_x - B^x_y - E^x_z \right) \\
\dot B^y = \frac12 \left( E^z_x + B^y_x - B^x_y - E^x_z \right) \\
\dot B^z = \frac12 \left( -E^y_x + B^z_x - E^x_y - B^x_z \right)
\end{cases} \\
x = 1 ~ \begin{cases}
\dot E^y = \frac12 \left( -E^y_x - B^z_x + E^x_y + B^x_z \right) \\
\dot E^z = \frac12 \left( -E^z_x + B^y_x - B^x_y + E^x_z \right) \\
\dot B^y = \frac12 \left( E^z_x - B^y_x + B^x_y - E^x_z \right) \\
\dot B^z = \frac12 \left( -E^y_x - B^z_x + E^x_y + B^x_z \right)
\end{cases}
\end{gather}
Applying a similar treatment to $y$ and $z$ bounds,
\begin{gather}
y = 0 ~ \begin{cases}
\dot E^x = \frac12 \left( E^x_y + B^z_y - E^y_x - B^y_z \right) \\
\dot E^z = \frac12 \left( E^z_y - B^x_y + B^y_x - E^y_z \right) \\
\dot B^x = \frac12 \left( -E^z_y + B^x_y - B^y_x + E^y_z \right) \\
\dot B^z = \frac12 \left( E^x_y + B^z_y - E^y_x - B^y_z \right)
\end{cases} \\
y = 1 ~ \begin{cases}
\dot E^x = \frac12 \left( -E^x_y + B^z_y + E^y_x - B^y_z \right) \\
\dot E^z = \frac12 \left( -E^z_y - B^x_y + B^y_x + E^y_z \right) \\
\dot B^x = \frac12 \left( -E^z_y - B^x_y + B^y_x + E^y_z \right) \\
\dot B^z = \frac12 \left( E^x_y - B^z_y - E^y_x + B^y_z \right)
\end{cases}
\end{gather}








\begin{gather}
z = 0 ~ \begin{cases}
\dot E^x = \frac12 \left( -(B^y+E^x)_z + E^z_x + B^z_y \right) \\
\dot E^y = \frac12 \left( (B^x-E^y)_z - B^z_x + E^z_y \right) \\
\dot B^x = -\dot E^y \\
\dot B^z = \dot E^x
\end{cases} \\
z = 1 ~ \begin{cases}
\dot E^x = \frac12 \left( (E^x - B^y)_z - E^z_x + B^z_y \right) \\
\dot E^y = \frac12 \left( (E^x + B^y)_z - B^z_x - E^z_y \right) \\
\dot B^x = \dot E^y \\
\dot B^y = -\dot E^x
\end{cases}
\end{gather}

To handle the edges, we must do three separate coordinate transformations: rotation about the $z$-axis (so up/down edges stay up/down), rotation about the $x$ axis (so left/right edges stay left/right), and rotation about the $y$ axis (so front/back edges stay front/back). Let's start with a rotation about the $z$ axis.

The following transformation will rotate the coordinate system $45^\circ$ about the $z$-axis counterclockwise, preserving the orthonormality of the coordinate vectors
\begin{gather}
x = \frac1{\sqrt 2} \left(r-s\right) \\
y = \frac1{\sqrt 2} \left(r+s\right) \\
z = t
\end{gather}
By preserving orthonormality, then in this coordinate system Maxwell's equations are still
\begin{gather}
\dot{\ura E} = \nabla \times \ura B \\
\dot{\ura B} = -\nabla \times \ura E
\end{gather}
This is important as we can now reuse the analysis above for the $(x,y,z)$ coordinates for these $(r,s,t)$ coordinates! In particular, this means that
\begin{gather}
r = 0 ~ \begin{cases} 
\dot E^s = \frac12 \left( -(E^s+B^t)_r + E^r_s + B^r_t \right) \\
\dot E^t = \frac12 \left( (B^s-E^t)_r - B^r_s + E^r_t \right) \\
\dot B^s = -\dot E^t \\
\dot B^t = \dot E^s
\end{cases} \\
r = \sqrt 2 ~ \begin{cases}
\dot E^s = \frac12 \left( (E^s + B^t)_r - E^r_s - B^r_t \right) \\
\dot E^t = \frac12 \left( (E^t + B^s)_r - B^r_s - E^r_t \right) \\
\dot B^s = \dot E^t \\
\dot B^t = \dot E^s
\end{cases}
\end{gather}
\begin{gather}
s = -\frac1{\sqrt2} ~ \begin{cases}
\dot E^r = \frac12 \left( (B^t - E^r)_s + E^s_r - B^s_t \right) \\
\dot E^t = \frac12 \left( -(B^r + E^t)_s + B^s_r + E^s_t \right) \\
\dot B^r = \dot E^t \\
\dot B^t = -\dot E^r
\end{cases} \\
s = \frac1{\sqrt2} ~ \begin{cases}
\dot E^r = \frac12 \left( (E^r + B^t)_s - E^s_r - B^s_t \right) \\
\dot E^t = \frac12 \left( (E^t - B^r)_s + B^s_r - E^s_t \right) \\
\dot B^r = -\dot E^t \\
\dot B^t = \dot B^r
\end{cases}
\end{gather}
Using that the set of points defined by $(0,0,z)$, $(1,0,z)$, $(0,1,z)$, and $(1,1,z)$ in $(x,y,z)$ correspond to $(0,0,t)$, $(\frac1{\sqrt2},-\frac1{\sqrt2},t)$, $(\frac1{\sqrt2},\frac1{\sqrt2},t)$, and $(\sqrt2,0,t)$, respectively. So we get that
\begin{gather}
(0,0,z) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + (B^x_z - B^z_x) + 3 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + (B^z_y - B^y_z) + 3 (B^z_x - B^x_z) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^y_z - E^z_y) + \sqrt2 (E^x_z - E^z_x) + 2 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^y_x - B^x_y) + (E^z_x - E^x_z) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^x_y - B^y_x) + (E^y_z - E^z_y) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^y_z - B^z_y) + \sqrt2 (B^x_z - B^z_x) + 2 (E^x_y - E^y_x) \right)
\end{cases} \\
(0,1,z) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_y - E^y_x) + (B^z_x - B^x_z) + 3 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^x_y - E^y_x) + (B^y_z - B^z_y) + 3 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + \sqrt2 (E^x_z - E^z_x) + 2 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^x_y - B^y_x) + (E^x_z - E^z_x) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^x_y - B^y_x) + (E^z_y - E^y_z) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^z_y - B^y_z) + \sqrt2 (B^x_z - B^z_x) + 2 (E^x_y - E^y_x) \right)
\end{cases}
\end{gather}
\begin{gather}
(1,0,z) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + (B^z_x - B^x_z) + 3 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + (B^y_z - B^z_y) + 3 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^y_z - E^z_y) + \sqrt2 (E^z_x - E^x_z) + 2 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^y_x - B^x_y) + (E^x_z - E^z_x) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^y_x - B^x_y) + (E^z_y - E^y_z) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^y_z - B^z_y) + \sqrt2 (B^z_x - B^x_z) + 2 (E^x_y - E^y_x) \right)
\end{cases} \\
(1,1,z) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_y - E^y_x) + (B^z_y - B^y_z) + 3 (B^x_z - B^z_x) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + (B^x_z - B^z_x) + 3 (B^z_y - B^y_z) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + \sqrt2 (E^z_x - E^x_z) + 2 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^x_y - B^y_x) + (E^z_x - E^x_z) + 3 ( E^y_z - E^z_y ) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^y_x - B^x_y) + (E^y_z - E^z_y) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^z_y - B^y_z) + \sqrt2 (B^z_x - B^x_z) + 2 (E^y_x - E^x_y) \right)
\end{cases}
\end{gather}

For the left/right edges, we'll make the transformation
\begin{gather}
x = t \\
y = \frac1{\sqrt2} (r-s) \\
z = \frac1{\sqrt2} (r+s)
\end{gather}
which will preserve
\begin{gather}
\dot{\ura E} = \nabla \times \ura B \\
\dot{\ura B} = -\nabla \times \ura E
\end{gather}
and this will continue into the exact same results as above, except with $(x,y,z)$ cycled to $(y,z,x)$. So our net result is
\begin{gather}
(x,0,0) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^z_x - E^x_z) + \sqrt2 (E^y_x - E^x_y) + 2 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + (B^y_x - B^x_y) + 3 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + (B^x_z - B^z_x) + 3 (B^x_y - B^y_x) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^z_x - B^x_z) + \sqrt2 (B^y_x - B^x_y) + 2 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^z_y - B^y_z) + (E^x_y - E^y_x) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^y_z - B^z_y) + (E^z_x - E^x_z) + 3 (E^x_y - E^y_x) \right)
\end{cases} \\
(x,0,1) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + \sqrt2 (E^y_x - E^x_y) + 2 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_z - E^z_y) + (B^x_y - B^y_x) + 3 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^y_z - E^z_y) + (B^z_x - B^x_z) + 3 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^x_z - B^z_x) + \sqrt2 (B^y_x - B^x_y) + 2 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^y_z - B^z_y) + (E^y_x - E^x_y) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^y_z - B^z_y) + (E^x_z - E^z_x) + 3 (E^x_y - E^y_x) \right)
\end{cases}
\end{gather}
\begin{gather}
(x,1,0) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^z_x - E^x_z) + \sqrt2 (E^x_y - E^y_x) + 2 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + (B^x_y - B^y_x) + 3 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + (B^z_x - B^x_z) + 3 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^z_x - B^x_z) + \sqrt2 (B^x_y - B^y_x) + 2 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^z_y - B^y_z) + (E^y_x - E^x_y) + 3 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^z_y - B^y_z) + (E^x_z - E^z_x) + 3 (E^x_y - E^y_x) \right)
\end{cases} \\
(x,1,1) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + \sqrt2 (E^x_y - E^y_x) + 2 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_z - E^z_y) + (B^x_z - B^z_x) + 3 (B^y_x - B^x_y) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_y - E^y_z) + (B^y_x - B^x_y) + 3 (B^x_z - B^z_x) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^x_z - B^z_x) + \sqrt2 (B^x_y - B^y_x) + 2 (E^z_y - E^y_z) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^y_z - B^z_y) + (E^x_y - E^y_x) + 3 ( E^z_x - E^x_z ) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^z_y - B^y_z) + (E^z_x - E^x_z) + 3 (E^x_y - E^y_x) \right)
\end{cases}
\end{gather}

Continuing this pattern for the front/back edges,
\begin{gather}
(0,y,0) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + (B^y_x - B^x_y) + 3 (B^y_z - B^z_y) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^x_y - E^y_x) + \sqrt2 (E^z_y - E^y_z) + 2 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + (B^z_y - B^y_z) + 3 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^z_x - B^x_z) + (E^x_y - E^y_x) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^x_y - B^y_x) + \sqrt2 (B^z_y - B^y_z) + 2 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^x_z - B^z_x) + (E^y_z - E^z_y) + 3 (E^x_y - E^y_x) \right)
\end{cases} \\
(1,y,0) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^z_x - E^x_z) + (B^x_y - B^y_x) + 3 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + \sqrt2 (E^z_y - E^y_z) + 2 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_x - E^x_z) + (B^y_z - B^z_y) + 3 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^z_x - B^x_z) + (E^y_x - E^x_y) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^y_x - B^x_y) + \sqrt2 (B^z_y - B^y_z) + 2 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^z_x - B^x_z) + (E^z_y - E^y_z) + 3 (E^x_y - E^y_x) \right)
\end{cases}
\end{gather}
\begin{gather}
(0,y,1) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + (B^x_y - B^y_x) + 3 (B^z_y - B^y_z) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^x_y - E^y_x) + \sqrt2 (E^y_z - E^z_y) + 2 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + (B^y_z - B^z_y) + 3 (B^y_x - B^x_y) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^x_z - B^z_x) + (E^y_x - E^x_y) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^x_y - B^y_x) + \sqrt2 (B^y_z - B^z_y) + 2 (E^z_x - E^x_z) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^x_z - B^z_x) + (E^z_y - E^y_z) + 3 (E^x_y - E^y_x) \right)
\end{cases} \\
(1,y,1) ~ \begin{cases}
\dot E^x = \frac14 \left( \sqrt2 (E^x_z - E^z_x) + (B^z_y - B^y_z) + 3 (B^y_x - B^x_y) \right) \\
\dot E^y = \frac14 \left( \sqrt2 (E^y_x - E^x_y) + \sqrt2 (E^y_z - E^z_y) + 2 (B^x_z - B^z_x) \right) \\
\dot E^z = \frac14 \left( \sqrt2 (E^z_x - E^x_z) + (B^y_x - B^x_y) + 3 (B^z_y - B^y_z) \right) \\
\dot B^x = \frac14 \left( \sqrt2 (B^x_z - B^z_x) + (E^x_y - E^y_x) + 3 (E^y_z - E^z_y) \right) \\
\dot B^y = \frac14 \left( \sqrt2 (B^y_x - B^x_y) + \sqrt2 (B^y_z - B^z_y) + 2 (E^x_z - E^z_x) \right) \\
\dot B^z = \frac14 \left( \sqrt2 (B^z_x - B^x_z) + (E^y_z - E^z_y) + 3 ( E^x_y - E^y_x ) \right)
\end{cases}
\end{gather}

\section{Constraint-free Maxwell's equations}

All previous cases were merely practice to find a way to enforce conservation of energy in numerical calculations of the constraint-free Maxwell's equations
\begin{gather}
\dot \phi = -k \phi + \nabla \cdot \ura E - 4 \pi \rho \\
\dot \psi = -k \psi + \nabla \cdot \ura B \\
\dot{\ura E} = K \ura E + \left(\nabla + \ura a\right) \times \ura B + \nabla \phi - 4\pi \ura J \\
\dot{\ura B} = K \ura B - \left(\nabla + \ura a\right) \times \ura E + \nabla \psi
\end{gather}
These will be significantly more complicated, being eight coupled equations instead of six. In matrix form,
%\begin{gather}
%\ura u = \left( \begin{array}{c} \phi \\ \psi \\ E^x \\ E^y \\ E^z \\ B^x \\ B^y \\ B^z \end{array} \right) \\
%\mathbb X = \left( \begin{array}{cccccccc} 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0



\end{document}